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Abstract 

We introduce a variational approach for the Quantum Inverse Scattering Method 
to exactly solve a class of Hamiltonians via Bethe ansatz methods. We undertake 
this in a manner which does not rely on any prior knowledge of integrability through 
the existence of a set of conserved operators. The procedure is conducted in the 
framework of Hamiltonians describing the crossover between the low-temperature 
phenomena of superconductivity, in the Bardeen-Cooper-Schrieffer (BCS) theory, 
and Bose-Einstein condensation (BEC). The Hamiltonians considered describe sys- 
tems with interacting Cooper pairs and a bosonic degree of freedom. We obtain 
general exact solvability requirements which include seven subcases which have pre- 
viously appeared in the literature. 

1 Introduction 

The Quantum Inverse Scattering Method (QISM) [1-3] is used for constructing quantum 
Hamiltonians with multiple conserved operators, and in turn obtaining their exact solu- 
tions by Bethe ansatz methods. A solution of the Yang-Baxter equation [4-6] may be 
used to construct a transfer matrix which generates the conserved operators of the Hamil- 
tonian. The Bethe ansatz used to obtain the exact solution can assume many forms. The 
original approach due to Bethe [7] is commonly termed the co-ordinate Bethe ansatz, 
whereas a more modern approach in the framework of the QISM (subject to a suitable 
reference state) is the algebraic Bethe ansatz [1-3]. However there are several variants 
between these two formats, which rely on functional relations and assumed analyticity 
properties of the eigenvalue spectrum to obtain the exact solution [8-12]. We note that 
in principle the implementation of the co-ordinate Bethe ansatz is not dependent on any 
prior knowledge of an associated solution of the Yang-Baxter equation, nor the conserved 
operators of the Hamiltonian that it generates. 

Progress in cold atom physics has yielded many studies into the nature of the crossover 
between the low-temperature phenomena of superconductivity from the Bardeen-Cooper- 
Schrieffer (BCS) theory and Bose-Einstein condensation (BEC) [13-16]. Early theoretical 
accounts emphasised the need to study Hamiltonians which explicitly incorporate coupling 
between Cooper pairs of atoms and bosonic molecular modes [17, 18]. The goal of the work 
presented here is to implement a strategy motivated by both the co-ordinate and algebraic 



1 



Bethe ansatz approaches to obtain a general class of exactly solvable Hamiltonians, with 
both a bosonic and Cooper pairing degrees of freedom, such that they model BCS-BEC 
crossover behaviour. As we will show, this approach gives a unified construction for classes 
of exactly solvable Hamiltonians with multiple free coupling parameters. It reproduces 
some seven exactly solvable subcases which have previously appeared in the literature 
(including some for which the bosonic degree of freedom is decoupled from the system, 
leaving a BCS Hamiltonian) [19-25]. The availability of such exact results makes amenable 
the computation of correlation functions [23,24,26] which potentially can be compared to 
experimental results. 

To highlight the mathematical motivation for our work, it is instructive to examine 
the evolution of integrable models of correlated electrons in one-dimension. Here there 
are many examples whereby models were first introduced and solved exactly via the co- 
ordinate Bethe ansatz approach, with the rederivation of the Hamiltonian through the 
Yang-Baxter equation coming as a later development. The first of these is the Hubbard 
model which was exactly solved by Lieb and Wu in 1968 [27]. It was not until 1986 that 
Shastry [28] first made connection to the Yang-Baxter equation for this model, with subse- 
quent understanding of the algebraic structure developing through the 1990s (e.g. [29,30]). 
The strong coupling limit of the Hubbard model leads to the t-J model, which was shown 
to be exactly solvable for particular choices of the coupling parameters through the works 
of Schlottmann [31] and Sarkar [32]. In this case it was not so long before the conserved 
operators were constructed through the Yang-Baxter equation [33,34]. Following on from 
this there was a flourish of activity throughout the 1990s in the study of exactly solvable 
models for correlated electrons via co-ordinate Bethe ansatz, which included the Bariev 
model [35], the g-deformation of the supersymmetric U model [36], and the Alcaraz-Bariev 
model which ultimately included both Hubbard and a solvable t-J model as particular 
subcases [37]. All three Hamiltonians were later shown to be derivable from a solution 
of the Yang-Baxter equation. Some works showing the connection to the Yang-Baxter 
equation, viz. [38, 39] for the Bariev model, and [40] for the g-deformed supersymmet- 
ric U model, were reasonably rapid developments. The algebraic formulations for the 
Bethe ansatz solution were given in [41-43]. Recently the origin of the solution of the 
Yang-Baxter equation for the Alcaraz-Bariev model has started to become clear [44]. 

There are two key differences between the construction of correlated electrons referred 
to above and our approach below. The first is that in the co-ordinate Bethe ansatz 
approach, the ansatz wave-functions for the correlated electron models are taken to be 
superpositions of plane waves. We will instead formulate the wave-functions as a factoris- 
able operator acting on a reference state, which follows the spirit of the algebraic Bethe 
ansatz [1-3] and relates back to Richardson's original calculations for pairing Hamiltoni- 
ans [19]. However we differ from the algebraic Bethe ansatz approach through a second 
key point in that we will work directly with a variational Hamiltonian. We follow the 
spirit of the co-ordinate Bethe ansatz approach [7,27,31,32,35-37] in that we do not re- 
sort to any knowledge of a transfer matrix or a set of conserved operators. By combining 
these two aspects of the co-ordinate and algebraic Bethe ansatz methods, we are able to 
formulate a powerful technique for constructing exactly solvable models in a very general 
fashion. 

The structure of the manuscript is as follows. We begin Section 2 by introducing a 
general Hamiltonian describing a reduced BCS model coupled to a bosonic degree of free- 
dom. In Section 2.1 we apply the Bethe ansatz to obtain constraints on the Hamiltonian's 
coupling parameters which are sufficient for exact solvability of the Hamiltonian. There 
are two cases which are dealt with separately in Section 2.1.1 and Section 2.1.2. In section 
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3 we discuss the connection between seven known exactly solvable subcases that are limits 
of the general model we consider. 

2 Variational Hamiltonian 

We consider the general family of pairing Hamiltonians coupled to a bosonic degree of 
freedom 

H = H -H 1 (1) 

where 

L 

H = aN + KiV 2 + f(zk)N k , (2) 
fe=i 

ILL 

Hi = (^Yg{z k )b bl + (5Y,^i^ b lh + ^Yg(z k )'^)blb s , (3) 

k=l k=l k,s 

for some complex-valued functions f(z) and g(z), and real-valued a, k, and a, which will 
be subject to certain solvability constraints yet to be determined. The overline notation 
denotes complex conjugation, which imposes that the above Hamiltonian is hermitian. 
The operators b k = c_ fe c|. and N k = b\b k for k > are hard-core Cooper pair creation and 
number operators, where c\. are fermion creation operators and Cooper pairs are assumed 
to consist of paired fermions of zero total momentum. Spin labels have been suppressed, 
so we do not imply anything about the spin properties (e.g., singlet, triplet) of the pairing. 
We also have a single bosonic mode with operators 6q and N . The particle operators 
satisfy the following commutation relations: 

[b j} b k ] =0 = [6t, 6+] Vj,k>0, 

M*]=0 Vj,k>0,jj:k, (4) 

MJ]=I, [h,bl] = I-2N k Vk>0. 

The Hamiltonian commutes with the total number operator N = No + ^2 k N k . Hamil- 
tonians in this family describe a system of bosonic molecules, condensed into a single 
bosonic degree of freedom, coupled to L Cooper pairs which are bosonic-like Cooper pairs 
that must observe an exclusion principle. The Hamiltonians consist of two parts. The 
diagonal part H Q , given in equation (2), describes the bosonic mode and allowed Cooper 
pair energy levels. The self-interaction term kNq describes a shift in the frequency of the 
bosonic mode as it is populated. The cross-interaction part Hi, given in equation (3), 
describes level dependent molecule-pair coupling and pair-pair couplings when either of 
j3 or a are non-zero. For the special case a = (3 = k = 0, which suppresses any action 
of the Hamiltonian on the bosonic part of the underlying Hilbert space, the Hamiltonian 
(1) reduces to a general form of the BCS Hamiltonian. It is for this reason we can refer 
to the general pairing Hamiltonian as a BEC-BCS crossover Hamiltonian. 

We do not expect that the Hamiltonian (1) is exactly solvable in general. However, we 
have found solvability conditions for various sub-classes of the Hamiltonian, in particular 
the cases i) k = of no self-interaction term and ii) a — of no BCS pair-pair scattering 
term. Taking appropriate limits, these two cases are shown to reduce consistently to the 
same Hamiltonian when k = = a. 
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2.1 Exact Solvability Constraints 

We would like to understand the extent to which an exact solution can be found for pairing 
Hamiltonians of the form given in equation (1) for the yet to be determined functions f(z) 
and g{z). To begin, motivated by the approach of Richardson [19], we assume the ansatz 1 

M 

i*>=n c fo)i°> 

3=1 

for the eigenstates of (1), where |0) denotes the vacuum state, 

L 

C(y) = 1 (y)bl + Y,h(y,z k )bl yeC 
k=i 

and h(y, z) is yet to be determined. We introduce the notation 



M 



M 



Using the identities in (4) the following commutation relations are found: 

KC(y)]=-r(v)i 

[b k ,C(y)] = h(y,z k )(I-2N k ), 
[N ,C(y)}= 1 (y)bl 
[N k ,C(y)] = h(y,z k )bl 
[NlC(y)] = 7 (i/)6j(I + 2JVo) I 



[H , C{y)\ = ai (y)bl + ^ f(z k )h(y, z k )b\ + ^(y)bl(I + 2N 
We note that for an operator O, 



k=l 



M 



o, n CM 



M / l-l 



M 



E(n^))[o,c^)] n 

3=1+1 



l=m \r=m 



Using this identity we can show 

M 

6o|*) 



3=1 

M 

E^-)i*i> 



|0> 



3=1 

M fj-1 

W = E (]l c (yr))[h,c( yj )} ' ' c( yi )\o) 



M 



3=1 \r=l / 1=3+1 

M M 

3=1 3,¥3 



(5) 



*For simplicity, throughout the manuscript we only consider eigenstates which do not involve blocked 
levels. For a review of the blocking effect we refer to [45]. There is no technical impediment to extend 
results to accommodate the general case with blocked levels, but we omit instances with blocked states 
for the sake of readability. 
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We may then calculate the following 

M ML M 

Af L Af L 

H 1 \*) = f3j2J2 rtvMzkHWj) + EEWMto. **)6jl*i> 
j=l fc=l j=l fc=l 

M L M L 

J,i^j k=l j=l k,s=l 

M L 

- o- Y 9{zk)g(z s )h{ yj , z s )h(y h z s )b\b\\m 3l ) . 

j,l+j k,s=l 

These equations give the action of the Hamiltonian H on the state In order to 

determine the exact solution we require that |^) be an eigenstate of H . Thus we look to 
solve the eigenvalue problem: 

(Ho-H,)^) =E\V) (6) 

for some scalar E. In general this will not be possible. The set of constraints required to 
find a solution to equation (6) are called solvability constraints and define the manifold in 
the coupling parameter space along which the Hamiltonian (1) sustains an exact solution. 

In order to find an exact solution we look to write the sums involving any of the vectors 
b\bl\tyji), b^bf.^ j{) , blbll^ji), b k \^fj), or bl\^ 3 ) as a combination of the vectors and 
1^}. As a first step we write the sum over b k bl\^ji) terms as a sum over a %b k \^ji) part 
and a b\.\^ 3 ) part. However, comparing the blb\\^/ji) and b b k \^ji) terms we see that the 
particular constraint chosen to reduce the blb\\^ji) terms must be compatible with one 
chosen to reduce the b Q b k \^ji) terms. For compatibility we introduce the constraint 



g(z s )h{y 3 ,z s )h(y h z s ) = k{y 3 ,yi)h{y h z s ) + k{y h y 3 )h(y 3 ,z s ), \/y h y 3 ,z s G C 
where k(yi,yj) is to be chosen later. This allows us to write 

ML ML 

j=l k=l j=l k=l 

ML ML 

~ 2 PJ2J2 k (w, yj)Kvj, + o- J2 Ya^aT^hiy,, z a )b\\*j) 

3,1+3 k=l j=l k,s 

M L 

~ 2a J2J29( z k) k (yi,y3) h (y3^ z s)blbl\^ 3 i). 

3,1+3 k,s 

We then utilise the definition (5) for |^) to express the sums in terms of the desired 
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vectors, 

ML ML M 

#il*> = PY, E 7(l/iM**)&Il*i> + EESMto. «*)6j|*i> - 2/3 E W)6j|*i> 
j=l fc=l j=l k=l 

M ML 

+ 2/3 ^ ^^^^(^^Sl^ + ^EE^^)^^'^)^!^) 

ML ML 

Consider the full expression 

M ML M 

H\0f) = (a + K)J2 7(Vi)&Jl*i> ~ EE^%- l*i> + 2/3 E J")&Sl*i> 

ML Mi ML 

+ E E f( z k) h (yj, - /3 EE^-M^li 1 ^) - a E E faM^iKvi, *.)&li*i> 

j=i fc=i i=i fc=i j=i fc,« 

ML ML 

+ 2cr E E^)Mm)^;>-2<7 ^j(«*)%,!ii)7(!li)i{6il*ji> 

M M 

-2 / 9 53*(y i) y,)7(w)656j|*ii> + « E rtvMvM\*fi)- 

At this point we observe the appearance of identical coefficients of the vectors ftty^lj) 
and a X^it=i 9 ( z k)b\.\ty j) for cr 7^ 7^ /3. In finding the exact solution we must choose 
constraints on the coefficients that are compatible with this. Thus, the next constraint 
must involve the terms that are not yet related. We seek to express the remaining sums 
over vectors b\.\^/j) as a combination of the vectors |\&), (3bl\^j) and o''^2^ =1 g(zk)bl\^j) 
such that the latter two have the same coefficients. We choose the following constraint to 
reduce the terms to the desired vectors 

f(z k )h(yj, z k ) = yjh(yj, z k ) + g(z k )r(yj) (7) 

h{yj, z k ) = _ y Vj/j, z k e C, 

where the r(yj) will be determined by compatibility requirements. With these choices we 
find 

ML ML ML 

^ £ ^)4l*i> = E E vMVi> **)&ll*i> + E E 0(**Ml/;)&Il*i> 

j=l fc=l j=l k=l j=l k=l 

MM ML 

j'=i j'=i i=i fc=i 
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Use of this relation leads to 

M M M L 

H\V) = J>|*) + $> + k- vMvM^j) ~ P **)&Sl*i> 

j=l j=l j=l k=l 

M M L 

+ 2/3 ^ /.•(//,, //,)/>„ x \' j ) + £5>(**)[Kl/i) - 07(l/;)]&Il*i> 

Af L ML 

~ a J2 ^9(zk)9&jh( yj , Ofe) + 2ff E E 9(zk)Hyj, yi)b\\^j) 

j=l k,s k=l 
M M 

M i 

fe=l 

At this stage we remark that the choice of constraint (7) sets the energy eigenvalues to a 

M 

standardised form E = ?/, . 

3=1 

Exact solvability is achieved by requiring that the coefficients of all terms other than 
|\&) cancel. To meet this requirement we choose (for g(z k ) ^ 0) 

L M 

(Vj - a- f^hiVj) + P E 9( z k)h(yj, z k ) = 2/3^2 HVji Hi) 

k=\ i+j 

L M 

PliVj) ~ r (Vj) + a Yl 9(zk)h(y j: z k ) = 2a KVv Vi) 
k=x m 

P (k(yj, yi)i{yi) + Hvh Vihivj)) = K i(yjh(yi) (8) 

& (k(yj, yi)i{yi) + k{y h yjh(yj)) = 0. (9) 

However, we note for k / / the constraints (8) and (9) are incompatible and we 
have at least two separate cases. We must now derive two sets of solvability conditions, 
one for each case, however, they should agree in the appropriate limit k = = a when 
they describe the same family of Hamiltonians. Some care will also be needed in the limit 
P ->■ 0. 

2.1.1 Case 1: No Pair-Pair Interaction (a = 0) 

The first case we look at is cr = 0. Here the solvability conditions reduce to 



g(zk)h(yj, z k )h(y h z k ) = k(y j} yi)h{y h z k ) + k{y h yj)h(y h z k ) (10) 

L M 

(yj - a- K)j(yj) + P Y 9{zk)h(yj, z k ) = 2(3 ^ KVi> Vi) ( n ) 

k=i ijtj 

Pl(yj) - r( yj ) = (12) 

P ( k (yj, yi)i{yi) + k{y h yj)i{y 3 )) = ^i(yjh(yi) (13) 
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ji^kj — ~T7 — \ 



Firstly, (12) is satisfied for non-trivial 7(2/^) only if 
Substituting (12) into (14) gives 

h(yj,Zk) 

and substituting this into (10) gives 

- T -^ (3g(z k )j(y :j ) (3g(z k )-f( yi ) _ M^Hm) ^g{z k ) 1 {y j ) 

9{ z k — T, \ 77 \ — K \yjiyi)—F7 \ r K{yi, Vj — 77 \ 

f( z k) - Vj f{zk> - yi f{zk) - yi f{zk) - yj 

Rearranging, 

^P(k(yj,yih(yi) + k(yi,yjh(yj)) -?-^ , , 0(k(yj,yih(yi)yj + k(yi,Vjh(Vj)yi 

J\ z k) 7 — r - 7 — ; 9{Zk)9{z k ) 



KvvVi) = it, \T(j/j)» 



i(yjh(yi) i(yjh(yd 

Since the right hand side does not depend on the parameter z k we require that, for some 
constants c\ and c 2 , 

P(Hyj,yih(yi)yj + KvuvdiivM = c ii(yjh(yi)> 
P(Hvvyih(yi) + Hyhyjhiyj)) = wivMvi), 

C2f(z k ) - (3 2 g(z k )g(z k ) = c x . 
However, compatibility with constraint (13) requires 

C 2 = K. 

The solution is thus 

ci - c 2 yi 

P(yj - yi) 

f(z k ) = eg 1 (fi 2 g{z k )g(z k ) + , 

C 2 = K, 

for some constant C\. 

Substituting this into the constraint (11) completes the compatibility of constraints 
yielding 

%•-(« + /€)+>—— = 2 > ( 15 ) 

C 2 = K. 

We refer to (15) as the Bethe ansatz equations for which the sub-family of Hamiltonians 
of (1) corresponding to a = are solvable. An equivalent expression is obtained by 
manipulating the terms, 

-L(a + K + c 2 (2M-2-L))-l + ]j2—/ 77 w 

2% 2 4 ^ yj /c 2 - f{z k )/c 2 

/ L 1 M \ M ( 16 ) 

+ [y - + y 2 = y 02 

2% V /(**) - j^vi-vi) i£ vo - vi 

C 2 = K. 
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In the special case of k = = a we find c\ = — (3 2 \g(z k )\ 2 = —G 2 , c<i = and the Bethe 
ansatz equations reduce to 



a yj/G 



1 2 - 1 



2G 2 A^f^/G-yi/G 2^y j /G-y l /G' 

We summarise by listing the constraining relations for exact solvability in the case 
a = here: 



/(z fc ) = « 1 /3 2 g(z k )g(z k ) + « x ci 

C 2 = «. 

for constants /3, and ci. These constraints define the manifold in the coupling parameters 
of the pairing Hamiltonian (1) for which it is exactly solvable with eigenstates of the form 
(5) in the case a = 0. In the above derivation we find that j(yj) is not fixed, however, it 
will be fixed by normalisation of the eigenstate. 

2.1.2 Case 2: No Self-Interaction Term (k = 0) 

Setting k = in the solvability conditions leads to the following set of constraints, 



g(zk)h(yj,z k )h(y h z k ) = k(yj,yi)h(yi,z k ) +k(y l ,y j )h(y j ,z k ) (17) 

L M 

(Vj - a h(Vj) + 9(zk)h(yj, z k ) = 2/3 k(y j: yi) (18) 
k=i i+j 

L M 

hiVj) - r(yj) + 9{zk)HVj, z k ) = 2a ^ k(y j: y t ) (19) 
fc=i i^j 

k(yj, yi)l{yi) + k{y h yjhiyj) = (20) 

Constraints (18) and (19) are compatible when 

Pr(yj) = [a{a - Vj ) + (3 2 }i{y 3 ) 
which along with condition (20) results in the constraint 

(3[a(a - Vj ) + f3 2 )k{ yj , yi )r{ Vl ) + 0[a{a - Vl ) + (3 2 ]k(y h yj )r( yj ) = 0. (22) 
and along with (21) results in 



Ph{y h z k ) = [a(a- yj ) + /3 2 ] 



2i#(^h(%0 
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Substituting (21) into (17) gives 



g{zk)r{yj) g{z k )r(yi) g(z k )r(yi) g(z k )r{y^ 



9i. z k)-FT-\ TT~\ = k (yj^yi)-FT~\ ^HVhVj 



3) 



f{z k ) - Vj f(z k ) - yi ' f(z k ) - yi ' f(z k ) - y 3 - 

which we rearrange to 

,/ ^(yj,yi) r (yi) + Kyhyj) r (yi) -t-t , ^ Hyj,yi) r (yi)yj + KvuVjYiyfiyi 

J\ z k) / — \ / \ 9{z k g{z k = — . 

r \yj) r {yi) r{y ] )r{yi) 

Since the right hand side of the equation does not depend on the parameter z k we must 
satisfy, for some constants c 2 and Ci, the following set of relations: 

HVj,Vi)r(yi)yj + k(yi,yj)r(yj)yi = cir{y j )r{yi), 
k(yj,yi)r(yi) + k{y h yj)r(yj) = c 2 r(y j )r(y l ) 1 
C2f(z k ) - g(z k )g(z k ) = c x . 

For compatibility of the first two equations we find the solution must be of the form 

yj yi 

However, for compatibility with constraint (22) we then require 

P[a(a - yj) + (3 2 ]^^r( yj )r( yi ) - P[a{a - W ) + /3 2 ] Cl ~ ^ r( yi )r( yj ) = 

Vj - yi yj - yi 

which we can rearrange to obtain 

(3[c 2 (aa + (3 2 ) - Cia}r(y j )r(y{) = 

and we must have 

c 2 (aa + (3 2 )f3 = ckj/3. (23) 
At this point we have the conditions 

cicr/3 = c 2 (aa + /3 2 )/3, 

Hvj'Vi) = - — — r (yj)i 

Vj - Vi 



ci = c 2 f(z k ) - g(z k )g(z k ), 
Prfy) = [a (a - yj) + /3 2 }l(yj), 
9{zk)r(yj) 



h(yj,z k 



For j3 7^ 0, these equations, along with equation (18), or equivalently (19), complete the 
compatibility of constraints yielding 

Vj - <* c 2 f(z k ) - Ci _ o Cl ~ C2Vl 

H*-y 3 ) + P 2 } f{z k )-y 3 f£ vs-vt ' 

Cl a = c 2 (aa + (3 2 ). (24) 



10 



For the case — 0, the constraint (24) is no longer necessary. In this instance the 
bosonic degree of freedom decouples and we may project onto a BCS system which will 
be discussed in the next subsection. 

An equivalent expression for the Bethe ansatz equations is obtained by manipulating 
the terms, 

a , A 1 



C { m-2-L))- — i— , + 



2 Vj \ (7 ( c i/ c 2 - Vj) ) 2a{cx/o2 - yj) A ^y j /c 2 - f(z k )/c 2 

/ L , M \ M 



c 2 (aa + (3 2 )(3 = C\G0. 



In the special case of a = = k we find c 2 = 0, (3 2 ci = —f3 2 g(zk)g(zk) = —G 2 and the 
Bethe ansatz equations reduce to 

a yj/G 1 



1 A 2 _ 1 

A 2^ t(y, \in -nJH ~ 9 



2G 2 i^ffayG-yj/G 2^ yj /G- m /G 

This result is in complete agreement with that of the first case, where we derived the 
Bethe ansatz equations for k^O and a = and then set k = at the end to obtain the 
limiting case. 



2.2 Recovering Known Exactly Solvable Subcases 

In the above subsections we determined manifolds in the coupling parameters of (1) for 
which an exact solution exists 2 . Taking appropriate limits of the general exactly solvable 
models yields eight subcases which we have presented in Figure 1. Seven of these subcases 
are known [19-25] 3 . 

The three tiers apparent in the graph in Figure 1 correspond to the number of free 
parameters in the solvable models which determine the coupling interaction strengths 
and the functional relationship between f(z) and g(z). The top tier consists of the most 
general exactly solvable models which were derived in sections 2.1.1 and 2.1.2. These 
models each have 3 free parameters up to an arbitrary energy rescaling. Models in the 
middle tier have 2 free parameters and models in lowest tier are described by only 1 free 
parameter. 



2 Whilc wc have tried to be general in deriving this manifold, it is potentially possible to relax some 
of the assumptions made on the functions introduced in Section 2.1. This will be considered further in 
future work. 

3 In some of the subcases a change of variable is required to bring the Hamiltonian into the form cited. 
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BCS/BEC [24] 



BCS/BEC [23] 



BCS [f] 



a=n(L-2M) 




BCS/BEC [21] 



Dicke [20] 




BCS [22] BCS [25] BCS [19] 



Figure 1: Hierarchy of known exactly solvable models that are limiting cases of our 
solutions for a = and k = 0. Key: [24] is the 2-channel p + ip-w&ve BCS coupled to 
bosonic molecular pairs with self-interaction term and no pair-pair interaction; [23] is the 
p + ip-w&ve BCS coupled to bosonic molecular pairs with no self-interaction term; [21] is a 
BCS system coupled to a bosonic mode which was derived by use of Sklyanin's boundary 
QISM; [20] is equivalent to the Dicke model as studied by Gaudin; [22] is the p + ip- 
wave BCS model; [25] is a BCS model with an energy cut-off used for the study of heavy 
nuclei; [19] is the original reduced s-wave BCS case due to Richardson; [I] is a general 
exactly solvable BCS model which it appears has not previously been explicitly identified 
in the literature. 

For the a = case the Hamiltonian is of the form 

L L 

H = aN + kN 2 + f(*k)N k (<7(*k)&o&l + K^bth 



k=l 

and is exactly solvable with eigenstates 

M 



k=l 



\ f{z k ) - Vj 



bl |0> 



when 



(3 2 \g(z k )\ 2 + d = c 2 f(zk), c 2 = n 



M 



2 E 



c\ - Kyi 



* Vj-Vi 



(25) 



Equation (25) affects the pairing symmetry in the BCS part of the model. The 3 inde- 
pendent parameters in this model are j3, k, and c\. 
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For the k = case the Hamiltonian is of the form 



L L L 

k=l k=l k,s 

and is exactly solvable with eigenstates 
when 

c 2 (cxa + /3 2 ) = cicr, 

\g(z k )\ 2 + Cl = c 2 f(z k ), (26) 

c ifa- - g) | c 2 f(z k ) - ci _ 2 Ci - c 2 yi 
c 2 - c i%' ~ f( z k) ~ Vj ^ Vj - Vi ' 

Immediately we notice the similarity with the a = case in that the equation (26) affects 
the pairing symmetry in the BCS part of the model. 

In the limit of j3 = 0, the k = model reduces to a model with no coupling between 
the bosonic degree of freedom and the Cooper pairs. Here we can simply project out the 
bosonic degree fo freedom from the Hilbert space of states and consider the Hamiltonian 
as BCS type only. The BCS model in this case has neither p+ip-wave or s-wave symmetry, 
but a generalisation of both. The Hamiltonian is of the form 

L L 
k=l k,s 

and is exactly solvable with eigenstates 4 




!*>HlllJ^f^|lo> 



where 



1 L 



\g{z k )\ 2 + c 1 = c 2 f{z k ), (28) 

L ,/ , M 

C2j{Zk) ~ Ci _ Ci ~ C 2 Vl , 2Q , 

Here C\ = results in p + ip-wave symmetry, while c 2 = results in s-wave symmetry 5 . 

We take this moment to reflect on some of the historical passage of events in integrable 
pairing Hamiltonians. Although the construction of exact eigenstates in the s-wave model 
traces back to Richardson's work of 1963, it was only in 1997 that the conserved operators 



4r 



These eigenstates are obtained by renormalising the eigenstate (5) to accomodate for the factor 

i 

singlet type in the s-wave case. 



lX=i "fiVj) which approaches zero in the limit /3 — > 0. 

5 Strictly speaking we also have to impose that the pairing is triplet type in the p + ip-wave case and 
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were constructed [46]. It was soon realised that these operators could be reproduced fol- 
lowing the approach of Gaudin [20] by considering the rational class of Gaudin magnets. 
By extending this to the trigonometric and hyperbolic classes more general integrable 
Hamiltonians were constructed in 2001 [47,48]. It turns out that the conserved operators 
in the hyperbolic case are those of the p + ip-w&ve Hamiltonian, a model which came to 
prominence through the work of Read and Green in 2000 [49]. However it was not imme- 
diately apparent that the hyperbolic class of conserved operators could be combined to 
produce the p+ip-wave Hamiltonian, and it was only in 2010 that the explicit relationship 
was made [26]. 

It is somewhat surprising that the above results show that the rational and hyperbolic 
cases can be seen as two limits of a more general pairing Hamiltonian (27) with the 
constraint (28) and the Bethe ansatz equations (29). The number of free parameters is 
the same in both limiting Hamiltonians. This is in some contrast to the case of the XXX 
spin chain, which is the rational limit of the XX Z spin chain obtained by a one-variable 
reduction in the coupling co-efficient of the S^S^ +l interactions. 

3 Conclusion 

We introduced a variational approach for the Quantum Inverse Scattering Method to 
exactly solve a class of Hamiltonians via Bethe ansatz methods. The procedure was 
conducted in the framework of variational Hamiltonians describing BCS-BEC crossover 
physics through interacting Cooper pairs and a bosonic degree of freedom. We obtained 
general exact solvability requirements which included seven subcases which have previ- 
ously appeared in the literature. An initial question to consider is whether the general 
forms of exactly solvable Hamiltionians, in Subsection 2.1.1 and 2.1.2 respectively, do 
admit a set of conserved operators. We have found that such a set is indeed constructible, 
either through the Gaudin algebra approach of [50] (by relaxing constraints imposed in 
Subsection 3.1 of their work), or by a suitable adaptation of the classical Yang-Baxter 
equation approach of [51]. However we emphasize that our motivation was to undertake 
calculations to obtain exactly solvable Hamiltonians in a manner which did not rely on 
any prior knowledge of integrability through the existence of a set of conserved operators. 

Finally, we remark that there remains scope to further extend this approach to a more 
general level than that which we have considered here. For example the exactly solvable 
Russian doll BCS model [52, 53] does not fit into the above scheme, both in that the 
pair-pair interaction is not factorisable, and the wave-function ansatz is of a different 
type. Furthermore pairing Hamiltonians such as those in [47,48,51] contain interaction 
terms which are not present in our starting Hamiltonian (1), and similar extended exactly 
solvable models can also be constructed within the above approach. 
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